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Chapter 1 


Introduction 


Among various types of domain decomposition methods, the work of M.E. 
Rose et al. on initial-boundary value problems by finite volume methods is gaining 
well-deserved attention. This approach offers a fresh look at century-old problems 
arising in the mechanics of continua. As can be seen from the literature review at 
the end of this report, the growing body of that effort is not yet available to wider 
readership. Cited work is primarily theoretical, whereas here some two-dimensional 
problems arising in solid and fluid dynamics are solved in conjuncture with Rose s 
formulation and the presentation of the method itself is reduced to the necessary 
minimum. 

All the presented applications involve the numerical treatment of harmonic and 
bi-harmonic operators. Therefore, some simple facts about the Poisson equation 
have to be reviewed first. 

For a plethora of problems, the following equation is created: 

j>gradw*dS = F ( 1 - 1 ) 

It is customary to apply divergence the theorem and rewrite Eq. 1.1 in the form: 

J\7 2 wdV=F (1.2) 
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For two-dimensional rectangular grids this procedure leads to a standard five-point 
stencil representation: 

lj 2u>j j -f- + Ky( w i y j-1 “h fhj (1*3) 

where /c x = and K y = 

The right-hand side term fij is usually interpretted as a source (per unit volume 
or area) and is given at the center of each cell. 

However, this entire procedure can be easily modified for a new, specific defini- 
tion of the normal derivative. 

For reasons to be explained shortly, consider another function 4> — $( x i y)- ft is 
assumed that the computational domain T> is divided into I * J non-overlapping, 
rectangular cells and the symbol <f>ij refers to the value of the function (f> at the 
center of ( ij ) -th cell. Whenever this argument of the function is displaced by half- 
width of cell ‘left’ ( ‘right’ ) the same function will be referred to as u,_j (u,+j). 
Similar notation is used when ‘upward’ or ‘downward’ displacement is considered. 
A simple diagram for the above notation is presented in Fig. 1.. 


u i,j+ 


u i-,j $i,j u i+,j 


u i,J- 


Fig.l. T>ij cell. 

Half-step normal derivatives are given by: 

( 1 . 4 ) 

( 1 . 5 ) 
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v h = 

V xij = - U i-j) 
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V 3U = /\y^ U ‘ ,1+ ~ 

*W = 


( 1 . 6 ) 

(1.7) 


It is assumed that cells T>ij form a non-overlapping cover of the entire compu- 
tational domain T> and it will be hereon required that functions u and derivatives 
v be continuous across centers of faces of neighboring cells. It was shown in [1] 
that values <f>ij play only auxiliary role, whereas the requirement of the continu- 
ity of normal derivatives furnishes additional equations necessary to complete the 
formulation. This requirement can be written in the form: 


^xi j 

— v xi+l,j 

y*j 

= v yi,j+i 


( 1 . 8 ) 

(1.9) 


Substituting u for w , and yet unspecified, function -W for F, the entire analysis 
is repeated for half-step derivatives 1.4-1. 7, for each cell T>ij. This time however, 
direct summation replaces the use of the divergence theorem. As a result the 
following system of equations is obtained: 

7x( u i-,j — 2 + u,+j) + 7 y( u i,j- — 2 $ ij + u i,j+) + w ij — (1-10) 


where: j x = 2k x and -y y = 2/c y . 

This should be complemented by equations stemming from the continuity of 
normal derivatives. 

The scheme introduced above is an example of a new, weak compact scheme, 
in terminology introduced by Rose , [1]. ‘New’ refers to an improved modeling of 
normal derivatives; the half-step formulas that were absent from [4], whereas ‘weak’ 
refers to the fact that only derivatives normal to faces of cells are involved. No com- 
prehensive study of the practical applications was presented to date. Some pilot 
results are given in [6]. Numerous problems arise when applying compact scheme 
in higher dimensions and/or systems of equations. Among them, the difficulty of 
matching this ‘stingy’ formulation with an equally efficient numerical procedure. In 
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the context of standard engineering problems of continuum mechanics, three differ- 
ent types of solutions were attempted. The first one, based on Paceman-Rachford 
alternative direction implicit (ADI) scheme [2], revealed serious stability problem. 
Second attempt involved the application of Kaczmarz algorithm [3] for the entire 
system of equations (including the ADI) and was hampered by slow covergence and 
low accuracy, although several potential advantages became apparent. Finally, a 
compromized approach based on the Gauss algorithm was tested and proved to be 
a viable alternative. 



Chapter 2 


Coupled Systems of Equations 


For the square domain ~D, we will consider a coupled system of time-dependent 
equations given by: 


I* = 

V 2 $ = -q>. 

We will impose the boundary conditions in the form: 

$(x, 0) = <£(:r, 1) = $(0, y ) = $(l,y) = 0 

and 


( 2 . 1 ) 

( 2 . 2 ) 


(2.3) 


&*<*■ °> = £«<*' » = = = ° W 

Initial conditions are given only for the ^-function by: 

$(0 ,x,y) = $o(z,y) (2.5) 

It is important here, that the boundary conditions are imposed only on the 
function $, therefore the method of solution considered further in this report, has 
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to take into the account the transfer of these boundary conditions. Considerations 
based on Green’s theorem lead to a simple resolution of this problem. 

Combining Eq’s 2.1 and 2.2 with 


— ^V 2 $)dV = * ndS 


the following relation is obtained: 

d 


J$~V<N+ J{V 2 + $f)(N = 

* ndS 


(2,6) 

(2.7) 


Furthermore, introducing the notation: 
v = V4> 
u = V<f> 

t>„ = * n and 

u n = V$ * n , respectively, and applying another form of the Green’s theorem: 


one arrives at: 


I V$* vrw + 


I $v 2 rw 


/ 


$vr * ndS 


(2.8) 


— j u 2 cW + j (* 2 + *f)<N = (2.9) 

j)[$(v n + — U„) - U n ^]dS. 

This relation, Eq. 2.9 above, constitutes the energy principle for the entire sys- 
tem. For the type of boundary conditions expressed by Eq. 2. 3-2. 4, the uniqueness 
of a solution for a steady case can be easily proven. 

Consider the forcing function, f, to be equal to zero everywhere in the com- 
putational domain, T). Then the function $ becomes a solution to the Laplace 
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equation. Moreover, the boundary conditions are uniformly equal zero and, as a 
result $ = 0, everywhere in T>. For the time-dependent problem, the analysis of 
Eq. 2.9 leads to the conclusion that the system has a contracting property, that is, 
solution converges asymptoticaly to zero, everywhere in T>. 

It is worth mentioning that the boundary conditions used here are related to 
the dissipative boundary conditions in the sense of Kreiss [8]. In fact, the left-hand 
side of Eq 2.9 needs to be non-positive for the assertion about the uniqueness of 
the solution to hold, just like in Kreiss conditions, except that his statements were 
applied to a single function relations. This observation will be very useful for some 
practical applications considered further in this work. 

Discrete Energy Principle 

In the context of compact schemes, the discrete equivalent of the energy prin- 
ciple was first proven for the diffusion equation by Rose [9] and [10]. Derivation, 
below, constitutes its extensions necessary for cases including coupled systems of 
diffusion and Poisson equations. 

The notation introduced for gradients will be used in rewriting the system of 
Eq’s 2. 1-2.2. in the form: 


v$ 

= u 

(2.10) 


= V 

(2.11) 

V * u + 

= 0 

(2.12) 

— T 

dt 

= V * v — /. 

(2.13) 


It is easier to follow the derivation when u and v are interpreted through their 
respective components, namely, u = (p, q) and v = (r,s). 

It will be assumed that the Cartesian grid is uniform with steps Ax, and Ay in x- 
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and y-directions respectively. In the same way, A t is a step in time. 

Following [ 1 ], half-step formulas for derivatives are introduced, capitalized Greek 
letters refer to values of functions at centers of elements: 


ptj = 

Ax^ i+ *’ j 

- *«) 

(2.14) 

PTj = 

&*«- 

ti-y) 

(2.15) 


Ay^ iJ+ i 

- *«) 

(2.16) 

9« = 

=<*- 

hj-0 

(2.17) 

r+ = 
•J 


- *y) 

(2.18) 

r ij ~ 

h*'”~ 


(2.19) 

$+. — 
S *J 


- *w) 

( 2 . 20 ) 

S ij = 



( 2 . 21 ) 


Standard differencing and averaging operators ( in time and in space ) will be 
used extensively: 

f*T 7 = ^(7i+l +7,-i) ( 2 - 22 ) 

6x7 = (7i+$-7i-p ( 2 - 23 ) 

Similarity, for y-direction and time operations. In order to make the notation 
more transparent, “tilde” ( like in 7 ) and “t” subscript ( like in 7 % ) will be used 
for time averaging and time differencing, alternatively. 

Applying the // and 8 operators to the half-step formulas, Eq 2.14-2.21, as it 
was done in [ 1 ], important relations for finite difference operators are obtained: 


A xfi x p — 8 x (f> 


(2.24) 
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Ayn y q = 6y<}> 

(2.25) 

A xp x r — 8 x xj) 

(2.26) 

Ay^yS = Sylf) 

(2.27) 

Ax 

$ - p x (f> - 8 x p 

(2.28) 

Cn 

< h 
1 

II 

& 

(2.29) 

/\x 

$ = - 4 S x r 

(2.30) 

A y 

* = PyV - ~f 6 v s - 

(2.31) 


Indices were ignored in the above formulas. The very nature of differencing and 
averaging, presented here, produces additional identities, regardless which of the 
variables, x, y or t is involved: 

Kip) = Mt )Kp) + 6(j)Kp) ( 2 - 32 ) 

p) = p(~f)p(p) + ^(7)W- ( 2 - 33 ) 

The first of the identities, Eq. 2.32, constitutes the fulcrum of the proofs con- 
tained in [10], because of the simple fact that = ^( 7 )^( 7 )* This was 

followed by a suitable defintions for both terms on the right-hand side. 

As it was already mentioned before, the continuity of normal derivatives will be 
required for any two adjacent cells, e.g. pfj — P7+ij- This fact is essential in the 
formulation of boundary conditions corresponding to Eq. 2.4, but it also manifests 
itself in rather non-intuitive spatial relations, specific to the compact scheme. We 
will re-interpret the model system, Eqs. 2.1 and 2.2, using the discretization al- 
ready introduced. 


Equations: 


Ax 6xP,j + A y 6y9ij 


At 1 






(2.34) 


(2.35) 
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are consistent with the system 2.10-2.13. In fact, the introduction of Eqs. 2.34 and 
2.35 defines the way in which spatial and time relations will be handled. More- 
over, these equations, Eq’s 2.24-2.31, assumptions about the continuity of normal 
derivatives and relations 2.32-2.33 are sufficient to prove the discrete energy princi- 
ple. Only a sketch of this proof will be presented here, because necessary algebraic 
manipulations are tedious but fairly elementary. 

Consider the expression “ ”. Using the identities 2.24-2.31 in Eq. 

2.35 it can be proven that: 

+ ^^3/(^0 — Px(pr) — Py(q$) ~~ *$/• (2.36) 

There is, however, an additional relation obtained through the time-differencing 
of Eq. 2.34. Applying similar algebraic manipulations to this relation, one obtains: 

Wft) - Px[\{p 2 )t] - Py[\{q 2 )t] (2-37) 

Both relations, Eq 2.36 and Eq. 2.37, are valid for each cell in V. In order 
to finish the proof, terms // x (pf) and p y (qs) have to be eliminated. This can be 
accomplished by a formal manipulation of Eq. 2.34. First the relation: 

-^Wp) -^Sytyq) = -* 2 + fi x (pr) + p.y(qs) (2.38) 

is obtained, and a simple observation can be made. In spite of Eq. 3.33, Eq. 2.38 
can be rewritten for time-averaged quantities, because it was originally based on 
linear, spacial relations, that could have been time-averaged beforehand. 

Recombining Eq. 2.36 and Eq. 2.37 with Eq. 2.38 ( for time-averaged quanti- 
ties ), the formal equivalent of the energy principle per cell is obtained: 
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M(i/> 2 )<] + ^[(js 2 )<] + * 2 + */ = 

- ^»W«) + ^ W«‘) 


(2.39) 


In order to produce a more transparent form, Eq. 2.39 has to be written for 
each cell of the domain and then the summation can be be performed. Some of the 
terms on the right-hand side cancel due to the requirement of continuity of normal 
derivatives, as well as the continuities of functions themselves, which is enforced 
implicitly through the notation already introduced. In the most consise form and 
neglecting the forcing function, the final result is given by: 


^j«.IM| S +ll*lf= (240) 

E»=i + Pt) - fpp] Ay - 

Ei=i,...,J,j=i[<£(r + p t ) - 4 >p]Ay + 

Zi=l,j=i,...,jW + Qt) ~ $q\ Ax - 

S,=ij=i,...,j[^(r + q t ) - rpq] Ax. 

The analogy with the Eq. 2.9 is complete in the sense that integral norms 
are replaced by the discrete-summation ones. Moreover, the homogenous bound- 
ary conditions for the function $ and its normal derivatives are enforceable in a 
straightforward way. Therefore, the discrete solution is bound by the boundary 
and initial conditions and should converge well. In fact, second order accuracy 
is expected for the function, its derivatives and its Laplacian. The generalized 
dissipative boundary conditions expressed in Eq. 2.9 and Eq. 2.40, need not to be 
exactly homogenous. The proper statement should be limited to the requirement 
of non-positivness of the right-hand sides of these equations. This fact is particu- 
larily tempting, because the Eq. 2.1 is a diffusion equation by itself, whereas the 



13 


result of Eq. 2.40 holds regardless of what kind of boundary conditions are imposed 
on function 3'. This is true, regardless whether the specific implementation of a 
numerical scheme, at some solution step does or does not impose some additional 
numerical boundary conditions on this function. Furthermore, there should be a 
plethora of cases where the diffusion equation, by itself, should have a unique and 
well-behaved numerical solution for some regions in X>, thus forcing the uniqueness 
of the other, coupled function <3>. 

In relaxation of the requirements stated in Eq. 2.4, two points of view need to 
be taken into the account: 

• First, as it was proven earlier by Rose [ 9 ] ,[10] , for the diffusion equation, 
similar discrete energy principle can be derived. Therefore, the additional 
boundary conditions on T in some region, say Xfy,, is correct, provided that 
on the circumference of this domain &D y,, the dissipative conditions for the 
diffusion equation are met. Moreover, for a coupled system of equations, the 
dissipative nature of the boundary conditions for a system may be a result 
of a coupling and a particular shape of the boundary. 

• Second, for a number of important engineering problems, dissipative nature 
of the boundary conditions may sometimes be proven using fairly elementary 
physical arguments. 

Example 1. Static load on a rectangular plate. 

In Chapter 4, Eqs. 2.1 and 2.4 are used to describe deflection of a thin, square, 
homogenous plate under a static load / . The function $ is interpretted as a 
displacement. Homogenous boundary conditions represent the solid support, and 
the zeroing of normal derivatives is usually referred to as “clamping”. 

If no “clamp” is required on one of plate’s edges, say y — 0 or the ” bottom” 
edge, then in the vicinity, 

d 2 $ 

— vl/ = 

dy 2 


(2.41) 
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and 


_ 

dy 


(2.42) 


have the same sign, except for some very complicated distributions of the load. This 
is due to the fact that the u n term represents the slope while the other term, — T 


describes the attenuation of this slope away from the edge. Since the displacement 


of all edges is still assumed to be zero, boundary conditions for the entire system 
are dissipative. In this case, the imposition of “undamped” boundary condition for 
a specific edge is necessary, because the term — is senstitve to the orientation 


of the normal vector. 


Example 2. Driven cavity flow. 

In Chapter 5, modified Eq 2.1 and 2.3 are used to describe a flow in a square 
cavity. Functions and $ represent vorticity and stream function, respectively. 
Homogenous boundary conditions imposed on the stream function correspond to 
no-penetration requirement appropriate for solid boundaries and/or streamlines, 
whereas normal derivatives of the same function, represent components of velocity, 
tangent to the boundaries, which are zero, due to the no-slip requirement. It is as- 
sumed that three of the walls are stationary and the motion of the fluid is imposed 
by a steady movement of the third ’’wall” in its own plane. In agreement with the 
formulas: 





and 


(2.43) 



(2.44) 


value of u n represents the velocity of the “plate” driving the motion. At the same 
time, 4/ is the vorticity at a center of a cell adjacent to the moving boundary. It 
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is a matter of a simple verification to realize that u n and 4/ have always the same 
sign. Since homogenous boundary conditions for are still required, boundary 
conditions are dissipative. 

Example 4. Mixing motion driven by an agitator. 

Another example of fluid motion is of concern here, because it represents yet an- 
other way of imposing of boundary conditions on the system discussed in this 
report. It is assumed that on a circumference of a square domain P, no-slip and 
no-penetration boundary conditions are imposed. At the center of the domain, at 
centers of some centrally-located cells, a sinusoidal variation of vorticity is required. 
This emulates a motion of a “square” agitator/impeller of a mixer. 

The dissipative nature of boundary conditions is probably a result of the separa- 
tion of energy principles for the vorticity and stream functions along some unspeci- 
fied flow line. This can also be discussed using straightforward physical arguments. 
Consider the tern: 

ipUn = VK - v y n x + v x n y ). (2.45) 

In agreement with already introduced formulas, we will sum this term along the 
circumference of a square region V € surrounding all the cells representing the ag- 
itator. When vorticity 4/ has a positive value, then there will always be a region 
in which vorticity has a postitive value, due to imposed continuity of all functions. 
This holds, for grids which are fine enough and provided that it is the agitator that 
drives the fluid motion. Summation corresponding to the term : 

J ']>(- v y n x + v x n y )dS (2.46) 

reveals that it is always negative . The other part of the integral depends on the 
proximity of the solid walls and the strength of the mixer. That is, in some specific 
cases, its absolute value may be smaller than the other part. This itself suggests 
that under some circumstances, the solution to the flow ’’inside” the agitator may 
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be unique. As a result, the function $ which is the stream function has a unique 
solution in the region X\, where no additional conditions were ever imposed on this 
function. 

In order to prove the existence of the solution in the domain Dy = D — T> e , we 
will consider a streamline, surrounding the agitating cells, but still internal to the 
region T> s . Since no mass addition is allowed, such a streamline exists, therefore 
boundary conditions are dissipative along this line. Moreover, they are such along 
the outer circumference of the domain Dy- 



Chapter 3 

Note on Crank-Nicholson and ADI schemes 


As it was pointed in [1] and in [9], the compact scheme for the diffusion equa- 
tion in Cartesian elements can be reduced to a Crank-Nicholson or ADI method 
involving only the center-cell values. 

Using the earlier-required continuities of normal derivatives, which apply only at 
the inner boundaries between elements, the potential forms are obtained. 


Consider the following condition: 

r r+ij = r i ( 3 * 1 ) 

which in conjunction with Eq. 2.18 and Eq. 2.19 implies: 

= + (3-2) 

This procedure may be applied to any of the variables involved in the system , as 
well as to the derivatives, c.f. Eq’s 2.14-2.21. For example: 

r* = (3 3) 

t -*<-'.>> (3 ' 4) 

For a coupled system of equations, Eqs. 2. 1-2.2, boundary conditions are imposed 
on the function <£, therefore for the Poisson subsystem the original formulation is 
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kept. However, the diffusion subsystem may be cast in a lot simpler form. Referring 
specificaly to the Eq. 2.35, it will be noted that the time-stepping operation involves 
only center-cell values of function 4/, provided that the right-hand side could be 
rewritten in a similar fashion. This is accomplished using identities 3.3 and 3.4. 
Further, we will prove that the entire analysis following the Eq 2.40 still holds. 


• First, we will state that: 


Ax Sxnj + A y 6ySij ~ 


(3.5) 




holds for any “inner” ij-th Cartesian element. This relation has the form 
of the Crank-Nicholson scheme, except that it requires the imposition of 
boundary conditions at centers of cells, adjacent to the boundary of the com- 
putational domain. For the coupled system of equations, however, this issue 
may be resolved by an inspection of Eq. 2.34. It will be noted, that for 
the Laplace subsystem, all normal derivatives were retained, therefore for all 
boundary cells, conditions 2.3 and 2.4 are enforceable. That way, conditions 
needed for the values of at the centers of these cells are obtained. 


• Second, except for the boundary cells, relations 2.39 hold, together with the 
continuity of functions and normal derivatives. Therefore, Eq. 2.40 is valid, 
provided that the original formulation is retained for all cells adjacent to the 
boundary. 

Let’s, accept such a scenario; that is, assume that Eqs. 2.33 and 2.34 are 
valid there. Now, the off-boundary values of normal derivatives can be re- 
placed using the formulas similar to Eqs. 3.3 and 3.4, while the off-boundary 
values of $ can be replaced using Eqs. 3.1 and 3.2. Furthermore, as it was 
already noted before, for each of the boundary cells an equation describing 
the center-cell value was already obtained, therefore by subtraction, addi- 
tional equations for values of normal derivatives at all boundaries would be 
produced. However, using the half-step formulas Eqs. 2.18-2.21 and obtained 
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center-cell values of *1', “edge” values of this function would be immediately 
produced. Finally, except for the corner cells, where some caution has to 
be exerted, relations for the boundary values of function 'F and its normal 
derivatives could be solved a posteriori , permitting one to handle the diffu- 
sion subsystem as a Crank-Nicholson scheme. 

For corner cells, the conditions stated by Eq’s 2.3 and 2.4 will be strictly 
enforced. This allows for a straightforward analysis of edge values of function 

V. 

As it was indicated in [2] the alternative direction implicit method [ ADI ] of- 
fers a very interesting alternative to the Crank-Nicholson scheme. It has been 
applied to Stokes-type flows before, in [2] and [7]. In [6] the scheme was proven 
to be adaptable to the compact scheme formulation. For a two-dimensional 
formulation, ADI may be best described as a half-time-step method with x- 
and y- differentiations frozen alternatively, resulting in equations being tridi- 
agonal. 

Unfortunately, for coupled systems of equations, this scheme reveals high sen- 
sitivity to boundary conditions. In fact, a very complicated backtracking is 
usually necessary in order to assure the numerical stability. It would be very 
hard to review how this numerical instability is offset by different numerical 
gadgets. We will limit ourselves to only one general statement. Refering to 
Eq. 2.39, it is safe to say that in handling x- and y- differentiation at differ- 
ent time levels, there is simply no method to guarantee that the boundary 
conditions remain dissipative at the same time level. 



Chapter 4 


Applications in Solid Mechanics 


Numerical solutions of two systems of equations were analyzed for efficiency 
and accuracy. The first system, Eq’s 4. 1-4.2 below, describes a steady deflec- 
tion of a thin, square, clamped plate, under a distributed load. The second 
system, Eq’s 4.6-4. 7, is an extension of the former one, where additional 
time dependence is included. In a rigorous sense, the second system does 
not describe the time evolution of the solution to a thin plate problem, it 
however serves as a model for other problems in continuum mechanics, like 
fluid dynamics, chemical flows and heat transfer. It is particularly suited for 
two-dimensional Stokes flows. Some examples of practical interest in fluid 
dynamics will be considered in the next Chapter. 

Steady Problem 

The steady problem, may be solved directly using the transfer of boundary 
conditions described already. All the existence and uniqueness conditions 
hold just like in the unsteady case. However, the system itself reduces to a 
pair of coupled Poisson equations and is not very interesting. Therefore as far 
as numerics is concerned, all the steady cases will be treated as asymptotic 
(in time) cases of the unsteady formulation. 
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Consider three functions u = u(x, y), w(x, y) and /i(z,y) fullfilling the fol- 
lowing system of equations: 


V 2 w = fx 

(4-1) 

3 

1 

II 

> 

(4.2) 


for 0 < x < 1 and 0 < y < 1 

This system is usually rewritten in terms of a single fourth-order equation: 

V 4 u = f\ (4.3) 

where typical boundary conditions are: 

u(x, 0) = u(x, 1) - u(0, y ) = u(l, y) = 0 (4.4) 

and 

^ U (x, 0) = 1) = ») = 0 < 4 - 5 ) 

In general, system of Eq 4. 1-4.2 with boundary conditions of Eq 4. 4-4. 5 may 
be cast in terms of two eliptic equations with Dirichlet-type boundary con- 
ditions. Then, the coupling is accomplished through the forcing term in the 
second, and through the boundary conditions for the first equation. Specific 
formulation of the boundary conditions depends on the numerical strategy 
to be used. 

Solution to the steady problem, may be viewed as an asymptotic solution to 
the unsteady problem, presented next. 
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Unsteady Problem 

Time-dependence is introduced ad hoc in the following way: 

ttw = V 2 w - f 2 (4.6) 

dt 

V 2 u = -w (4.7) 

Together with boundary conditions Eq 2.4 and initial condition: 

u(0,x,y) = u 0 (x,y ) (4.8) 

Boundary conditions for Eq, 4.6 are specific to the strategy and are discussed 
as a part of the discretization process. Herein we will only point out to the 
discussion following the Eq. 3.5. It was indicated there that the conditions 
of the type: 

P\j = o (4.9) 

are sufficient to produce the boundary relations for Eq. 4.8. 

Discretization Process 

The square [0, 1] <8> [0, 1] is divided into 1 2 square cells, where indices i,j refer 
to cell’s centers. Time, t, is indexed with the superscript, n, and the following 
notation is used: 


t 1 

= o, 

(4.10) 

t U+ 2 

= t n + 0.5 * At, 


t N 

Tjjiax 



= w(t n ,Xi,yj), 

(4.11) 

ft 

= u(t ,r,, t/j), 

(4.12) 

n 

= u(t n ,Xi + 0.5Ax, yj — 0.5 Ay). 

(4.13) 
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Paceman-Rachford-type alternative direction implicit (ADI) scheme [2] for 
Eq. 3.6. takes the form: 


w 


+V“?J-l - 2 <i + <*l) - 0 . 5 Af/ 2j j 

n +h 

n . . * 


n+i n \ / n +2 o i n +2 \ 

ij “ w iJ — ^ Wjj + ij) 

(4.14) 


n-H 2 


Wi 


:r-u,:p= a , («#( - 2 ^ + 

'” +S (4.15) 


n-H 


n- 1- J 


+A s «+\ - 2 <+' + wifi,) - 0 . 5 A «/S’, 


where coefficients A x and A y are given by: 


0 . 5 A* 
Ax * Ax 
0 . 5 A* 
Ay * Ay 


(4.16) 

(4.17) 


Boundary conditions for Eq’s 4.14 and 4.15 are not known apriori therefore 
they have to be calculated at each step from the solution to the scheme for 
Eq. 4.7. As already mentioned, the Poisson problem Eq. 4.7 is cast in terms 
of a Compact Finite Volume Method [1]. In the notation used previously, 
this scheme can be rewritten in the following way: 


< + 7 x «-j ~ 2 + «&,■) + 7„(«S- - 2C, + «S+) 

<t>T- i,i-2<y + C> 

*S-i- 2 «u- + *3 


0, (4.18) 
0, (4.19) 
0. (4.20) 


For square cells, Ax = Ay, it is prudent to take 7^ = 7 y = 1 rather than 
7x = and 7 y = This procedure is equivalent to the scaling of 

the w function and it improves the conditionning of the entire system. Now 
on, this rescaling will be implicitly assumed. 
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Boundary conditions for the Eq. 4.18, 4.19 and 4.20 are generated using 
Eq.4.4 and this subsystem may be viewed as a Dirichlet problem for the 
Poisson equation. In order to account for the boundary conditions expressed 
by Eq. 4.5, for cells adjacent to boundary of the computational domain one 
would welcome relations of the type: 


„ 771 7771 

u l-,j ~ <Pi,j 

= 0 

(4.21) 

„,77i l m 

u I+,j - 9lj 

= 0 

(4.22) 

<1- = Cl 

= 0 

(4.23) 

C/+ = hi 

= 0 

(4.24) 


New boundary conditions for the Eq. 4.14 and Eq.4.15. follow directly from 
the Eq. 4.21-4.24. and the definition of the compact scheme. Refering to the 
Eq. 4.2. specifically, for the lower ( y = 0 ) boundary, the following relations 
are obtained: 


C = “7x(C,i + u i+,i) ~ 7j/(«m+)- ( 4 - 25 ) 

Similar relations hold for all the cells adjacent to the boundaries of the com- 
putational domain and need not to be repeated here. Eqs. 4.5, although 
similar in form, are not substituted for any of Eqs. 4.18. They, together with 
the Eq. 4.14. and Eq. 4.15. constitute a Dirichlet problem for the diffusion 
equation. The entire system is now specified for a set of 21 * (21 +1) variables, 
per each half-time step m = 1, 1 + 2, 2 + \ ... 

ADI Routine 

The very first attempt to solve the problem was essentialy based on methods 
presented by Peyret and Taylor [2]. Tridiagonal solver for Eq’s 4.14-4.15 was 
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implemented. Similarity, Eqs. 4.18-4.20 were cast in time-dependent form and 
a specialized tridiagonal solver, accounting for the continuity relations, Eqs. 
4.19-4.20, was developed. Boundary conditions for the system were produced 
at each time half-step. Moreover, since Eq. 3.25 does not contain relaxation 
parameter, various methods of such inclusion were tested. The system was, 
nevertheless persistently unstable. 

Compact scheme does not benefit from explicit Taylor expansions, there- 
fore higher-order boundary conditions become rather convoluted. Except for 
methods of Forouk and Fusegi [7], one is still faced with the necessity of back- 
tracking in time in order to correct the mid-step boundary values. Rather 
than using these complicated methods of obtaining proper boundary condi- 
tions, a global approach based on successive overrelaxation routine (SOR) 
was tested. 


Kaczmarz Routine 

The system of algebraic equation, Eqs. 4.14-4.15 and Eqs. 4.18-4.20 together 
with the suitable boundary and initial conditions was solved using Kaczmarz 
Successive Over relaxation Scheme [3]. A similar approach was used for an 
earlier version of a Compact Scheme by Gatski, Grosch and Rose ( GGR 
code, [4]). As noted in [4], Kaczmarz algorithm is very easy to implement, 
because equations may be programmed in groups and, therefore amenable to 
the out-of-core computations one might add. For a particular system consid- 
ered herein, the A matrix is not only sparse but it has almost cellular form. 
It is quite easy to pivot just few rows in order to obtain a system which 
can be vectorized and synchronized (parallelism) to a very large extent. The 
most simple, although somewhat redundant example of such a procedure can 
be given by the following. By subtracting Eq. 4.25 from a suitably indexed 
Eq. 4.18, and using the fact that function u has zero values on the bound- 
ary, and carrying similar procedures for all cells adjacent to the boundary 
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of the computational domain, complicated boundary conditions expressed by 
equations similar to Eq. 4.25 may be replaced by <f>ij = 0. A dearth of similar 
strategies is possible, provided that available computer has at least four pro- 
cessors and a typically slow convergence of SOR algorithm does not hamper 
the computational process too much. Results presented in this subsection 
and parts of the next one were obtained on a 64- bit Convex Cl 20 computer 
working in a sequential mode. 

For each half-time step, all dependent variables can be renumerated and 
absorbed into vector X, with J — 21 * (21 + 1) entries. System can be 
symbolically written as: 


Ai* X = b t (4.26) 

where Ai is a row-vector corresponding to an i-th equation of the algebraic 
system. We will assume that the L 2 norm of Ai is always 1. 

Define the residue to be: resi = Ai * X — bi 

and introduce X° as the initial guess for the vector of unknowns, then the 
Kaczmarz procedure consists of the following: 

Step 1. X? +1 = X? - aA,res ?; for : i = 1, ... J — l,n = 0, 1,2, . . . 

Step 2. X? +1 = X] 

The procedure is repeated until the absolute value of residue res achieves 
required minimum. The value of the relaxation parameter a should vary 
between 1 and 2. 


Results I 

For all the testing some analytic functions had to be selected. They were 
constructed in the following way: The x- and y-dependence as well as time- 
dependence were introduced through X = x 2 *(1 — x) 2 , Y — y 2 *(1 — y ) 2 and 



27 


T = 0.5 * (1 — t) * (2 — t) and the function u was introduced by u = T * X *Y. 
Functions w and // ( ft ), needed for calculations and error norm evaluations, 
were obtained analytically by substitutions into Eq. 4.2 and Eq. 4.1 ( or 
Eq.4.6), respectively. Shapes of functions u, w and // are drawn in Fig.2. 


f function 
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w function 



function w 
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u function 



displacement u 
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Fig. 2 Functions u , w and /; . 
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Steady Problem 

The solution process of the steady problem involves the ADI scheme, therefore 
the ‘ t ’ variable had to be interpreted as ‘ficticious’ time and T = 1. In this 
case, the actual solution becomes an asymptotic (in time) solution to the 
unsteady problem and the initial guess becomes merely the initial condition. 
For all runs with ‘ficticious’ time, u 0 = 0 was selected, however in order to 
speed up the process a simple tridiagonal solver for Eqs. 4.14. and 4.15. was 
introduced ahead of the Kaczmarz routine. This solver was used only once 
for the very first half-time step. 

In order to establish the second-order accuracy of the method, maximum 
error and / 2 error norms for u and (unsealed) w functions were computed 
for various grid coarsness. Using the values of functions u and w obtained 
analitically following the norms were computed: 


I U max — vnax | u anal u comp | 


IV max — max W ana l lV com p | 


II U II 2 “ 2/*(2J+l)(£ ( w *-+,J-+ana/ u »-+,j- 

2 1 

II ^ || 2 = j(S (^ i } janal ^tjeomp) ) 


+ com P f) i 


Table l.a Error norms for u and w functions 


maximum residue = 10 5 


no.of cells || u Umax II u II 2 II W Wmax II W 11 2 
3 2 0.00160326 0.00171 0.02515275 0.01977 

6 2 0.00036485 0.00060 0.01268210 0.00610 

12 2 0.00010021 0.00021 0.00489993 0.00152 
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Table l.b Error norms for u and w functions 


maximum residue = 10 7 


no.of cells 
3 2 
6 2 
12 2 


u 


II « II; 


w 


w 


0.00160131 0.00171 0.02518925 0.01979 
0.00036470 0.00060 0.01264418 0.00608 
0.00009965 0.00021 0.00485032 0.00155 


As seen from the Table 1, above, the system gives a second-order accurate 
solution to the biharmonic equation. The accuracy of the auxiliary function 
w is lagging because of rescaling used to improve the conditioning of the 
entire system. 

The relaxation parameter a should be optimal, in the sense that it should re- 
duce to the minimum the number of iterations required to satisfy the residue 
criterion Experiments showed that the optimal value of the relaxation param- 
eter a varies strongly with the number of cells in the computational domain. 


Table 2. No of Cells versus optimal a 


3 2 

1.54 

5 2 

1.75 

6 2 

1.84 

10 2 

1.90 

12 2 

1.91 

20 2 

1.82 


Even for an optimal a, the number of iterations needed to achieve the re- 
quired minimum varied strongly with the number of cells and the value of 
the threshhold. 
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Table 3. No of Cells versus the number of iterations 
resi = 10 -5 ,re52 = 10~ 6 ,re.S3 = 10 -8 

no.of cells resi res 2 rcs$ 

5 2 62 72 123 

10 2 1016 1201 ... 

20 2 35628 

Results presented in Table 3. show that the number of necessary iterations 
grows rapidly with the number of cells. Moreover, it is inversely proportional to 
the degree of tolerance imposed on the residue of the system. 

Unsteady Problem 

Because of the rescaling of the w -function, solutions to the unsteady problem 
become very expensive; extremly low tolerance on residue has to be required in 
order to prevent errors from becoming much larger than the actual of values of the 
solution. Nevertheless, a steady problem can still be handled using a supercomputer 
with an advantage of significantly lower memory rquirements. 

Gauss Routine 

The idea of pivoting of the A matrix, introduced previously was brought to a 
natural conclusion by replacing the SOR algorithm with a simple Gauss transfor- 
mation with pivoting. The simple IBM matinv subroutine was adopted. The 
ADI scheme was also abandoned in favor of Crank-Nicholson scheme, for Eq. 3.6. 
For completeness it is given below: 


A, (<C?; - 2<+> + wgfj) 
+ A, «+!, -2<+‘ + W ?+i) 
+ a* he,., - 2"", + u, r+ij) 




+ A, -2< + <*,)- Af 


(4.27) 



34 


Introduction of the Gauss transformation allowed for a more efficient study 
of properties of solutions. Results presented in Table 1. were recalculated using 
matinv routine and are presented below: 

Table 4. Error norms for u and w functions 

steady case 


no. of cells 
3 2 
6 2 
12 2 


II ti IU, II U || 2 
0.00154321 0.00167 
0.00036470 0.00060 
0.00009965 0.00006 


II W Wmax II W 112 
0.02623457 0.02002 

0.01264307 0.00608 
0.00484978 0.00155 


The corresponding unsteady results are presented in Table 5. They were ob- 
tained using the same Gauss subroutine but on Cray Y-MP supercomputer. The 
upper entries correspond to the error values for time, t=1.0 and the lower entries 
correspond to the error values at t=2.0, cf. the definition of function T, in this 
subparagraph. 

Table 5. Error norms for u and w functions 

unsteady case 


no.of cells || u 

3 2 0.00000233 

0.00000195 
6 2 0.00000138 

0.00000125 
12 2 0.00000040 

0.00000036 


II w Umax 

II «> ||, 

0.00016804 

0.00006 

0.00014060 

0.00005 

0.00003912 

0.00002 

0.00003486 

0.00002 

0.00001284 

0.000007 

0.00001206 

0.000006 


Results presented in Table 5 show the ’’quality” of the approximation, since in 
the considered example values of functions u and w for t=l. 0,2.0 are zero. 




Chapter 5 

Applications in Fluid Dynamics 

In [2], the numerical treatment of a stream function - vorticity formulation 
for two-dimmensional Stokes flow was described in detail. Equations are usually 
written in the form: 



(5.1) 

0 = V 2 y + u 

(5.2) 

dx 

&X A 

oy 

(5.3) 

dx 

a » = ~Yx 

(5.4) 


In the equations above, u is vorticity, x - stream function, and a x and a y are 
components of velocity vector. Re, is the Reynolds number. 

Herein, an application of a weak compact scheme for handling of boundary 
conditions is considered. In general, whenever vorticity is involved in the problem, 
the strong version has to be implemented, because derivatives normal, as well as 
tangent to the faces of the cells are necessary. However, for a test case of wall- 
bound, vortical flows and a specific use of the scheme discussed here, it is still 
possible to employ only the less complicated, weak scheme and gain considerable 
savings in CPU time and memory. 
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The form of Eq. 5. 1-5.2 is essentially the same as Eq. 4.6-4. 7. Values of 
vorticity lj are evaluated at center-cell locations in the customary way. Second- 


order formulas for convection terms axe 

given by: 


(-> = 

Xi,j+ — Xi,j- 
Ay 

(5.5) 

(2*) = 

y dx’ij 

X*+,i Xi-j 
Ay 

(5.6) 

du) 

(—) = 
K dx } ij 

u i+ l,j ~ u i-lJ 
2Ax 

(5.7) 

(-) = 
1 dy ’a 

~ 1 

(5.8) 

2Ay 


A standard Newton-Raphson procedure was used to handle the nonlinear terms. 
In the final run the matinv subroutine was abandoned in favor of its fully vector- 
ized counterpart from Linpack [5]. Computations were carried out on Cray Y-MP 
supercomputer at North Carolina Supercomputing Center, Research Triangle Park, 
NC. 

Special care was taken to vectorize the code in order to make it efficient. Except 
for I/O operations and some complicated indexing routines, the code is thoroughly 
vectorized. Statistics showed that 60-87 percent of the time was spent in a subrou- 
tine rearanging Newton-Raphson iterations. 

Note on Upwind Differentiation 

In Chapter 3 the method of the transfer of boundary conditions was described in 
detail. The important facet of this methods was the elimination of “edge” values of 
the function 1 t. The entire procedure may be interpreted as a solution in the sense 
of staggered grids. Moreover, as it was shown in [6], such a formulation improves 
the accuracy. For fluid dynamics what is even more important, is the fact that 
there was no need to even consider values of vorticity at solid wall boundaries. In a 
sense, the center-cell value of vorticity is a resultant of mid-face values of velocity, 
c.f. Eq 2.34. This statement can, in principle be carried in full for the potential 
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form of the compact scheme . There axe, however some conceptual obstacles that 
need to be addressed before such a programme can be implemented. 

We will rewrite the system 5. 1-5.4 in a form: 


u = 



(5.9) 

V * u + \I> = 

0 


(5.10) 

v = 

eV^- 

- 

(5.11) 

II 

V * v 

-/ 

(5.12) 

a x — 

Uy 


(5.13) 

dy — 

—U x 


(5.14) 


System 5.10-5.15 is equivalent to Eqs. 5. 1-5.4, because Eqs. 5.9, 5.13 and 5.14 
guarantee that the divergence of the velocity is zero. Moreover, the energy princi- 
ple is readily available: 


ld_ 
2 dt 


J u 2 (N + e I * 2 dV + 
J $/dV - Jv(d* V$)dV = 

j,[$( — u„ + v n ) - €^u„]d5 


(5.15) 


Comparing the above equation with Eq. 2.9, we note that although there is a slight 
difference in the definition of the v vector, both equations have essentially the same 
form. There is an additional term in Eq. 5.15. Using Eq. 5.10 and 5.15 it can be 
proven that it is equal zero. Now, the entire analysis following the Eq. 2.9 may be 
easily repeated for the system 5.9-5.14. 

Surprisingly, this additional term which is handled so easily for the differential 
equation, causes rather severe problems in finite difference counterparts. This is 
due to the fact that for e = ^ >> 1 the system 5. 1-5.4 behaves quite like a 
hyperbolic one. In terms of numerical analysis, such a problem is resolved using 
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the “upwinding” or “upwind differentiation”. Methods for such a treatment were 
already developed in [10], but they were limited to a single diffusion equation. 
The essence of that modification consists of introduction of cell-Reynolds-number- 
dependent weights into formulas of the type 2.14-2.21, 

For a coupled system of equations, such as the one described in this report, “up- 
winding” has to be performed consistently for the stream function and the vorticity. 
Otherwise, due to the nature of the extra term in Eq. 5.16, no convergence could 
be guaranteed. In itself such a procedure is not very complicated. True difficul- 
ties arise when the upwind differentiation is employed in the study of wall-bound 
vortices. This is a topic for a separate study and only the approximate formulas 
5. 5-5. 9 will be used. Resulting numerical instabilities are not severe, they are easily 
pinpointed by the inspection of motion graphics, developed for the project. 

Results II 

Three practical problems presented below, have no known analytic solutions, 
therefore the analysis has to rely on a careful study of pictures, including the 
animations recorded on a video. Only results for 12x12 cell runs are presented. The 
choice of flows was not accidental. In order to keep things simple the advantage was 
taken from the fact that for the wall-bound vortices the value of stream function 
X on the circumference of the computational domain can be taken equal to zero 
without the loss of generality. 

Captured Wheel Flow. 

A two dimensional vortex is initially at a steady state defined by V$ = c * r . 
Vorticity has a constant, positive value. Instantaneously, impermeable walls are 
impressed upon the flow. Walls form a square concentric with the center of the 
vortex. In figures presented below, Reynolds number is defined as Re = — , 

where c is the proportionality constant of the vortex, L the length of the side of 
the bounding wall, v is the coefficient of kinematic viscosity. Boundary conditions 
used for this flow are standard for viscous flows. In the numerical context, they 
become the same as the ones used for the clamped plate, in the previous section. 
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Fig. 3 Captured Wheel Flow Re=10000 . 



Time step 0 



n 1 rt 1 2 


Captured Wheel Flow 
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Because of the continuous decay of the velocity field, no instabilities were de- 
tected. 

Driven Cavity Flow. 

A two-dimensional flow in a rectangular cavity, Fig. 4, is assumed to be viscous 
and initially at rest. The walls are impermeable. The velocity U is impressed upon 
the topmost layers of the fluid, instantenously. Reynolds number is based on the 
length of the gap and the magnitude of the impressed velocity. 




Time step 50 


ngapl 2 


driven Cavil y Flow 




Time step 100 


ngapl 2 


Driven Cavity Flow 
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Fig. 4 Driven Cavity Flow, velocity field at different times, Re — 10000 . 
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As in the previous case, setting the value of stream function x equal to zero at 
the circumference of the computational domain automatically enforces components 
of velocity normal to the wall to be zero. Components of velocity tangent to the 
boundary are defined as derivatives of x normal to the wall. Except the side open 
to the outside flow, all tangent components of velocity axe zero. As discussed in the 
paragraph on implementations of Kaczmarz algorithm, it is legitimate to impose 
conditions on normal derivatives of stream function, by manipulating the value of 
this function at centers of cells adjacent to the boundary. The result is written 
explicitly: 


Xil — — 0.5Ay * Ui (5.16) 

In this study only uniform values of velocity U are considered , but in general, 
any distribution of velocity is allowed. 

Corners of the computational domain, adjacent to the open flow, are singular 
points of the solution. Typically, various types of shifts are applied there, but 
for the presented scheme, experiments showed that almost any set of boundary 
conditions works well. It was decided to use the no slip boundary conditions in 
these two corners. 

Numerical instability is detectable after the influence of the upper plate reaches 
the right, bottom corner of the domain. 

Oscillatory Flow. 

Finally an oscillatory type of flow is presented. Here, a two-dimensional, vis- 
cous, unsteady flow is enforced by a periodic variation of the vorticity at the center 
of the domain. The intensity of the vorticity at the core is scaled with the frequency 
of oscillations / , and only one dimensionless parameter, Re — , appears in the 

equations. Boundary conditions at the wall are identical to the captured wheel 
flow. Initial conditions correspond to the stagnant flow. Fig. 5 represents various 
stages of the development of the oscillatory flow. 
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Fig. 5 Oscillatory Flow, velocity field at different time steps. 

Re = 10000. 
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In Fig. 6 the same flow, but for Re = 1000 is shown. The background, which is 
unfortunately given only in shades of gray, represents the intensity of the vorticity 
field. These pictures were transcribed from the color, motion graphics. The video 
was used to study in-phase, out-of-phase exctitation of the mixing motion. 
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Fig. 6 Oscillatory Flow, velocity field overlayed on the vorticity field, at differ- 
ent time steps. 

Re = 1000. 
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It was observed that the numerical instability becomes important well into the 
chaotic regime of the mixing flow. 

Discussion 

The aim of the presented study was to show the practical advantages of the 
compact scheme. All cases were limited to two-dimensional phenomena, although 
there were no such constraint in the method of solution. The code used to ob- 
tain presented results is somewhat limited by the specific formulation of the fluid 
dynamic problem. The method is fully explaineded in [1], where , the potential 
advantages of the compact scheme are discussed. However, before this poten- 
tial can be fully exploited, the computational efficiency has to be improved. The 
main issue is the characteristics of the matrix A of the system. There is a strong 
indication, [1], [6] and [10] that the compact scheme might benefit from the sym- 
metric, positive- definite matrices, although this avenue was not pursued in this 
report. However, the fact that it is now a coupled system of equations, defined on 
staggered grids plus some other aspects of the system of algebraic equations offer 
significant adavantages, that would have to explored sometime in the future. At 
the present time, however, the scheme offers a viable alternative to other, more 
popular methods. Finally, a chaotic flow was obtained using only 144 Cartesian 
elements, this best illustrates the potential of the scheme . 
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